#include"stdafx.h"
#include"hopnet.h"
#include"permut.h"

THopNet::THopNet(int New_n, int New_d, float New_eta) {
    Y = NULL;
    W = NULL;
    L = NULL;
    Permut = NULL;

    n = New_n;
    Setd(New_d);
    Seteta(New_eta);
    w = n*n;
    l = 0;
    if((Y = new float[n]) == NULL) {
        /* Error */
        n = w = l = 0;
        return;
        /* Error */
    } else {
        float* yi = Y;
        for(int i = 0; i < n; i++, yi++)
            *yi = -1.0;
    }
    if((W = new float*[n]) == NULL) {
        /* Error */
        n = w = l = 0;
        delete [] Y;
        return;
        /* Error */
    } else {
        float** wi = W;
        for(int i = 0; i < n; i++, wi++) {
            if((*wi = new float[n]) == NULL) {
                /* Error */
                n = w = l = 0;
                float** wj = W;
                for(int j = 0; j < i; j++, wj++)
                    delete [] *wj;
                delete [] W;
                delete [] Y;
                return;
                /* Error */
            } else {
                float* wij = *wi;
                for(int j = 0; j < n; j++, wij++)
                    *wij = 0.0;
            }
        }
    }
    if((L = new float*[l]) == NULL) {
        /* Error */
        n = w = l = 0;
        float** wj = W;
        for(int j = 0; j < n; j++, wj++)
            delete [] *wj;
        delete [] W;
        delete [] Y;
        return;
        /* Error */
    } else {
        float** lk = L;
        for(int k = 0; k < l; k++, lk++) {
            if ((*lk = new float[n]) == NULL) {
                /* Error */
                n = w = l = 0;
                float** lj = L;
                for(int j = 0; j < k; j++, lj++)
                    delete [] *lj;
                delete [] L;
                float** wj = W;
                for(int j = 0; j < n; j++, wj++)
                    delete [] *wj;
                delete [] W;
                delete [] Y;
                return;
                /* Error */
            } else {
                float* lki = *lk;
                for(int i = 0; i < n; i++, lki++)
                    *lki = 0.0;
            }
        }
    }
    Permut = new TPermut();
}

THopNet::~THopNet(void) {
    float** lk = L;
    for(int k = 0; k < l ; k++, lk++)
        delete [] *lk;
    delete [] L;
    float** wi = W;
    for(int i = 0; i < n ; i++, wi++)
        delete [] *wi;
    delete [] W;
    delete [] Y;
    l = w = n = 0;
    delete Permut;
}

int THopNet::Setd(int New_d) {
    if (New_d <= 0)
        New_d = 1;
    d = New_d;
    return d;
}

float THopNet::Seteta(float New_eta) {
    eta = New_eta;
    return eta;
}

int THopNet::Getn(void) {
    return n;
}

int THopNet::Getw(void) {
    return w;
}

int THopNet::Getl(void) {
    return l;
}

float* THopNet::GetY(void) {
    return Y;
}

float** THopNet::GetW(void) {
    return W;
}

float** THopNet::GetL(void) {
    return L;
}

int THopNet::Getd(void) {
    return d;
}

float THopNet::Geteta(void) {
    return eta;
}

int THopNet::ClearY(void) {
    #pragma omp parallel for
    for(int i = 0; i < n; i++)
        Y[i] = -1.0;
    return EXIT_SUCCESS;
}

int THopNet::ClearW(void) {
	#pragma omp parallel for
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n ; j++)
            W[i][j] = 0.0;
    return EXIT_SUCCESS;
}

int THopNet::ClearL(void) {
	#pragma omp parallel for
    for(int k = 0; k < l; k++)
        for(int i = 0; i < n ; i++)
            L[k][i] = 0.0;
    return EXIT_SUCCESS;
}

int THopNet::NegationY(void) {
    #pragma omp parallel for
    for(int i = 0; i < n; i++)
        Y[i] *= -1.0;
    return EXIT_SUCCESS;
}

int THopNet::LtoY(int pos) {
	#pragma omp parallel for
    for(int i = 0; i < n; i++)
        Y[i] = f(L[pos][i], d);
    return EXIT_SUCCESS;
}

int THopNet::RandY(int per) {
    int pn = per*n/100;
    int step = Permut->ReCreate(n,true);
    float* yi = &Y[step];
    for(int i = 0; i < pn; i++, yi+=step) {
        *yi = f((float)(RANDOM(100)-50)/50,d);
        step = Permut->GetStep();
    }
    return EXIT_SUCCESS;
}

int THopNet::RandW(void) {
	#pragma omp parallel for
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n ; j++)
            if(i != j)
                W[i][j] = W [j][i] = (float)(RANDOM(100)-50)/100;
            else
                W[i][j] = 0.0;
    return EXIT_SUCCESS;
}

int THopNet::RandL(void) {
	#pragma omp parallel for
    for(int k = 0; k < l; k++)
        for(int i = 0; i < n ; i++)
            L[k][i] = (float)(RANDOM(100)-50)/50;
    return EXIT_SUCCESS;
}

int THopNet::RandAllMatrix(void) {
    RandY(100);
    RandW();
    RandL();
    return EXIT_SUCCESS;
}

int THopNet::ReCreate(int New_n) {
	float** lk = L;
	for(int k = 0; k < l ; k++, lk++)
		delete [] *lk;
	delete [] L;
    float** wi = W;
    for(int i = 0; i < n ; i++, wi++)
        delete [] *wi;
    delete [] W;
    delete [] Y;
    n = New_n;
    w = n*n;
    l = 0;
    if((Y = new float[n]) == NULL) {
        /* Error */
        n = w = l = 0;
        return EXIT_FAILURE;
        /* Error */
    } else {
        float* yi = Y;
        for(int i = 0; i < n; i++, yi++)
            *yi = -1.0;
    }
    if((W = new float*[n]) == NULL) {
        /* Error */
        n = w = l = 0;
        delete [] Y;
        return EXIT_FAILURE;
        /* Error */
    } else {
        float** wi = W;
        for(int i = 0; i < n; i++, wi++) {
            if((*wi = new float[n]) == NULL) {
                /* Error */
                n = w = l = 0;
                float** wj = W;
                for(int j = 0; j < i; j++, wj++)
                    delete [] *wj;
                delete [] W;
                delete [] Y;
                return EXIT_FAILURE;
                /* Error */
            } else {
                float* wij = *wi;
                for(int j = 0; j < n; j++, wij++)
                    *wij = 0.0;
            }
        }
    }
    if((L = new float*[l]) == NULL) {
        /* Error */
        n = w = l = 0;
        float** wj = W;
        for(int j = 0; j < n; j++, wj++)
            delete [] *wj;
        delete [] W;
        delete [] Y;
        return EXIT_FAILURE;
        /* Error */
    } else {
        float** lk = L;
        for(int k = 0; k < l; k++, lk++) {
            if((*lk = new float[n]) == NULL) {
                /* Error */
                n = w = l = 0;
                float** lj = L;
                for(int j = 0; j < k; j++, lj++)
                    delete [] *lj;
                delete [] L;
                float** wj = W;
                for(int j = 0; j < n; j++, wj++)
                    delete [] *wj;
                delete [] W;
                delete [] Y;
                return EXIT_FAILURE;
                /* Error */
            } else {
                float* lki = *lk;
                for(int i = 0; i < n; i++, lki++)
                    *lki = 0.0;
            }
        }
    }
    return EXIT_SUCCESS;
}

int THopNet::AddL(void) {
    float** nL;
    if((nL = new float*[++l]) == NULL) {
        /* Error */
        l--;
        return EXIT_FAILURE;
        /* Error */
    } else {
        float** nlk = nL;
        float** lk = L;
        for(int k = 0; k < l; k++, nlk++) {
            if((*nlk = new float[n]) == NULL) {
                /* Error */
                l--;
                float** nlj = nL;
                for(int j = 0; j < k; j++)
                    delete [] *nlj;
                delete [] nL;
                return EXIT_FAILURE;
                /* Error */
            } else {
                float* nlki = *nlk;
                if(k<(l-1)) {
                    float* lki = *lk;
                    for(int i = 0; i < n; i++, nlki++, lki++)
                        *nlki = *lki;
                    delete [] *lk++;
                } else {
                    float* yi = Y;
                    for(int i = 0; i < n; i++, nlki++, yi++)
                        *nlki = *yi;
                }
            }
        }
        delete [] L;
        L = nL;
    }
    return EXIT_SUCCESS;
}

int THopNet::DelL(int dk) {
    if((l == 0) || (dk < 0) || (dk > l-1))
        return EXIT_FAILURE;
    float** nL;
    if((nL = new float*[l-1]) == NULL) {
        /* Error */
        return EXIT_FAILURE;
        /* Error */
    } else {
        float** nlk = nL;
        float** lk = L;
        for(int k = 0; k < l; k++, lk++) {
            if(k != dk) {
                if((*nlk = new float[n]) == NULL) {
                    /* Error */
                    float** nlj = nL;
                    for(int j = 0; j < k; j++)
                        delete [] *nlj;
                    delete [] nL;
                    return EXIT_FAILURE;
                    /* Error */
                } else {
                    float* nlki = *nlk;
                    float* lki = *lk;
                    for(int i = 0; i < n; i++, nlki++, lki++)
                        *nlki = *lki;
                    nlk++;
                }
            }
            delete [] *lk;
        }
        delete [] L;
        L = nL;
        l--;
    }
    return EXIT_SUCCESS;
}

int THopNet::PrintY(FILE* stream) {
    fprintf(stream,"N(%i)\n",n);
    for(int i = 0; i < n; i++)
            fprintf(stream,"%100.40f\n",Y[i]);
    return EXIT_SUCCESS;
}

int THopNet::PrintW(FILE* stream) {
    fprintf(stream,"W(%i)\n",n*n);
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n ; j++)
            fprintf(stream,"%100.40f\n",W[i][j]);
    }
    return EXIT_SUCCESS;
}

int THopNet::PrintL(FILE* stream) {
    fprintf(stream,"L(%i)\n",n*l);
    for(int k = 0; k < l; k++) {
        for(int i = 0; i < n ; i++)
                fprintf(stream,"%100.40f\n",L[k][i]);
    }
    return EXIT_SUCCESS;
}

int THopNet::PrintAllMatrix(FILE* stream) {
    PrintY(stream);
    PrintW(stream);
    PrintL(stream);
    return EXIT_SUCCESS;
}

int THopNet::ScanAllMatrix(FILE* stream) {
    float** lk = L;
    for(int k = 0; k < l ; k++, lk++)
        delete [] *lk;
    delete [] L;
    float** wi = W;
    for(int i = 0; i < n ; i++, wi++)
        delete [] *wi;
    delete [] W;
    delete [] Y;
    l = w = n = 0;

    int Scan_n = 0;
    int Scan_w = 0;
    int Scan_l = 0;

    fscanf_s(stream,"N(%i)\n",&Scan_n);
    n = Scan_n;
    if((Y = new float[n]) == NULL) {
        /* Error */
        n = 0;
        return EXIT_FAILURE;
        /* Error */
    } else {
        float* yi = Y;
        for(int i = 0; i < n; i++, yi++)
            fscanf_s(stream,"%f\n",&*yi);
    }
    fscanf_s(stream,"W(%i)\n",&Scan_w);
    if(Scan_w != n*n) {
        n = 0;
        delete [] Y;
        return EXIT_FAILURE;
    }
    w = Scan_w;
    if((W = new float*[n]) == NULL) {
        /* Error */
        n = w = l = 0;
        delete [] Y;
        return EXIT_FAILURE;
        /* Error */
    } else {
        float** wi = W;
        for(int i = 0; i < n; i++, wi++) {
            if((*wi = new float[n]) == NULL) {
                /* Error */
                n = w = l = 0;
                float** wj = W;
                for(int j = 0; j < i; j++, wj++)
                    delete [] *wj;
                delete [] W;
                delete [] Y;
                return EXIT_FAILURE;
                /* Error */
            } else {
                float* wij = *wi;
                for(int j = 0; j < n; j++, wij++)
                    fscanf_s(stream,"%f\n",&*wij);
            }
        }
    }
    fscanf_s(stream,"L(%i)\n",&Scan_l);
    if(Scan_l % n) {
        n = 0;
        delete [] Y;
        w = 0;
        for(int i = 0; i < n ; i++, wi++)
            delete [] *wi;
        delete [] W;
        return EXIT_FAILURE;
    }
    l = Scan_l/n;
    if((L = new float*[l]) == NULL) {
        /* Error */
        n = w = l = 0;
        float** wj = W;
        for(int j = 0; j < n; j++, wj++)
            delete [] *wj;
        delete [] W;
        delete [] Y;
        return EXIT_FAILURE;
        /* Error */
    } else {
        float** lk = L;
        for(int k = 0; k < l; k++, lk++) {
            if ((*lk = new float[n]) == NULL) {
                /* Error */
                n = w = l = 0;
                float** lj = L;
                for(int j = 0; j < k; j++, lj++)
                    delete [] *lj;
                delete [] L;
                float** wj = W;
                for(int j = 0; j < n; j++, wj++)
                    delete [] *wj;
                delete [] W;
                delete [] Y;
                return EXIT_FAILURE;
                /* Error */
            } else {
                float* lki = *lk;
                for(int i = 0; i < n; i++, lki++)
                    fscanf_s(stream,"%f\n",&*lki);
            }
        }
    }
    return EXIT_SUCCESS;
}

float THopNet::f(float s, int dl) {
    if(dl <= 0)
        dl = 1;
	double one_step = ((double)1.0)/((double)dl);
    if(s == 0.0)
		return (float)one_step;
    if(s > +1.0)
        return +1.0;
    if(s < -1.0)
        return -1.0;

	double step = ((double)s)/one_step;
	if (step > 0)
		return (float)(ceil(step) * one_step);
	if (step < 0)
		return (float)(floor(step) * one_step);
	return (float)(step * one_step);
}

float THopNet::MSE(float* N, float* nN) {
    float MSE = 0.0;
	#pragma omp parallel for reduction(+: MSE) shared(N, nN)
    for(int i = 0; i < n; i++) {
        MSE += (N[i] - nN[i])*(N[i] - nN[i]);
	}
    return sqrt(MSE/(float)n);
}

float THopNet::E(void) {
    float E = 0.0;
	float* pY = Y;
	float** pW = W;
	#pragma omp parallel for reduction(+: E) shared(pY, pW)
    for(int i = 0; i < n; i++) {
		for(int j = 0; j < n; j++) {
			E += (float)(W[i][j] * (Y[i] + 1.0) * (Y[j] + 1.0));
		}
	}
    E /= -8.0;
    return E;
}

float THopNet::CompNeuron(int pos) {
    float y = 0.0;
	#pragma omp parallel for reduction(+: y)
    for(int j = 0; j < n; j++) {
        y += W[j][pos] * Y[j];
    }
    y = f(y, d);
    return y;
}

float THopNet::MSECompNeuron(int pos) {
    float y = 0.0;
	#pragma omp parallel for reduction(+: y)
    for(int j = 0; j < n; j++) {
        y += W[j][pos] * Y[j];
    }
    y = f(y,d);
    float MSE = sqrt(((Y[pos] - y)*(Y[pos] - y))/n);
    Y[pos] = y;
    return MSE;
}

float THopNet::CompNeuron(float* pY, float** pW) {
    float y = 0.0;
    for(int j = 0; j < n; j++) {
        y += (*pW)[j] * Y[j];
    }
    y = f(y,d);
    return y;
}

float THopNet::MSECompNeuron(float* pY, float** pW) {
    float y = 0.0;
    for(int j = 0; j < n; j++) {
        y += (*pW)[j] * Y[j];
    }
    y = f(y,d);
    float MSE = sqrt(((*pY - y)*(*pY - y))/n);
    *pY = y;
    return MSE;
}

int THopNet::CompS() {
    float* nY;
    if((nY = new float[n]) == NULL) {
        /* Error */
        return EXIT_FAILURE;
        /* Error */
    } else {
		float* pY = Y;
		float** pW = W;
		#pragma omp parallel for shared(nY, pY, pW)
        for(int i = 0; i < n; i++) {
            nY[i] = CompNeuron(&pY[i],&pW[i]);
        }
		#pragma omp parallel for
        for(int i = 0; i < n; i++) {
            Y[i] = nY[i];
        }
        delete [] nY;
    }
    return EXIT_SUCCESS;
}

float THopNet::MSECompS() {
    float MSE = 0.0;
    float* nY;
    if((nY = new float[n]) == NULL) {
        /* Error */
        return EXIT_FAILURE;
        /* Error */
    } else {
		float* pY = Y;
		float** pW = W;
		#pragma omp parallel for reduction(+: MSE) shared(nY, pY, pW)
        for(int i = 0; i < n; i++) {
            nY[i] = CompNeuron(&pY[i],&pW[i]);
            MSE += (nY[i] - pY[i])*(nY[i] - pY[i]);
        }
		#pragma omp parallel for shared(nY, pY)
		for (int i = 0; i < n; i++) {
			pY[i] = nY[i];
		}
        delete [] nY;
    }
    return sqrt(MSE/(float)n);
}

int THopNet::CompA(void) {
	float* pY = Y;
	float** pW = W;
	#pragma omp parallel for shared(pY, pW)
    for(int i = 0; i < n; i++) {
        pY[i] = CompNeuron(&pY[i], &pW[i]);
    }
    return EXIT_SUCCESS;
}

float THopNet::MSECompA(void) {
    float MSE = 0.0;
    float y = 0.0;
	float* pY = Y;
	float** pW = W;
	#pragma omp parallel for reduction(+: MSE) shared(pY, pW)
    for(int i = 0; i < n; i++) {
        y = CompNeuron(&pY[i], &pW[i]);
        MSE += (y - pY[i])*(y - pY[i]);
        pY[i] = y;
    }
    return sqrt(MSE/(float)n);
}

int THopNet::CompR(void) {
	int i = Permut->ReCreate(n, true);
    for(int j = 0; j < n; j++) {
        Y[i] = CompNeuron(i);
        i += Permut->GetStep();
    }
    return EXIT_SUCCESS;
}

float THopNet::MSECompR(void) {
    float MSE = 0.0;
    float y;
    //int i = RANDOM(n);
	int i = Permut->ReCreate(n, true);	
    for(int j = 0; j < n; j++) {
        y = Y[i];
		Y[i] = CompNeuron(i);
        MSE += (y - Y[i])*(y - Y[i]);
        //i = RANDOM(n);
		i += Permut->GetStep();
    }
    return sqrt(MSE/(float)n);
}

int THopNet::CompMSES(void) {
    int c = 0, oc = 0, ocmax = n;
    float MSE =0.0f, oMSE = 0.1f, ooMSE = 0.2f, oooMSE = 0.3f;
    do {
        oooMSE = ooMSE;
        ooMSE = oMSE;
        oMSE = MSE;
        MSE = MSECompS();
        c++;
        if((ooMSE == MSE) && (oooMSE == oMSE)) {
            if(oc++ > ocmax)
                break;
        } else
            oc = 0;
    } while(MSE != 0.0);
    return c;
}

int THopNet::CompMSEA(void) {
    int c = 0;
    float MSE;
    do {
        MSE = MSECompA();
        c++;
    } while(MSE != 0.0);
    return c;
}

int THopNet::CompMSER(void) {
    int c = 0;
    float MSE;
    do {
        MSE = MSECompR();
        c++;
    } while(MSE != 0.0);
    return c;
}

int THopNet::LearnH(void) {
    if(l == 0)
        return 0;
    ClearW();	
    for(int k = 0; k < l; k++) {
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n; j++) {
                if(i != j)
                    W[i][j] += (float)((1.0/(float)n) * L[k][i] * L[k][j]);
            }
        }
    }
    return EXIT_SUCCESS;
}

float THopNet::LearnD(void) {
    if(l == 0)
        return 0.0;
    float MSE = 0.0;
    float* nL;
    if((nL = new float[n]) == NULL) {
        /* Error */
        return EXIT_FAILURE;
        /* Error */
    } else {
        float delta = 0.0;		
		float** pW = W;
		float** pL = L;
		float eta_n = eta/(float)n;
		int k = Permut->ReCreate(l, true);
        for(int d = 0; d < l; d++) {
			k += Permut->GetStep();
			#pragma omp parallel for shared(k, nL, pL, pW)
            for(int i = 0; i < n; i++) {
				nL[i] = 0.0;
				for(int j = 0; j < n; j++) {
					nL[i] += pW[i][j] * pL[k][j];
				}
				nL[i] = pL[k][i] - nL[i];
            }
			#pragma omp parallel for reduction(+: MSE) shared(eta_n, k, nL, pL, pW)
            for(int i = 0; i < n; i++) {				
                for(int j = 0; j < n; j++) {
                    if(i != j) {
                        pW[i][j] += eta_n * nL[i] * pL[k][j];
                        MSE += (eta_n * nL[i] * pL[k][j])*(eta_n * nL[i] * pL[k][j]);
                    }
                }
            }
        }
        MSE = sqrt(MSE/(float)(w*l));
        delete [] nL;
    }
    return MSE;
}

int THopNet::LearnDMSE(float MSEmin) {
    float MSE;
    int count = 0;
    ClearW();
    printf("Clear W matrix.\n");
    printf("Start LearnD function.\n");
    printf("#Step\tMSE\t\t(d=%i\teta=%1.10f\tlearnkind=d)\n",d,eta);
    do {
        MSE = LearnD();
        printf("%i\t%1.16f\n",count++,MSE);
    } while ((MSE > MSEmin) && (count < n));
    return EXIT_SUCCESS;
}

float THopNet::LearnMPR(void) {
    if(l == 0)
        return 0;
    float MSE = 0.0;
    float* oY;
    if((oY = new float[n]) == NULL) {
        /* Error */
        return EXIT_FAILURE;
        /* Error */
    } else {
        float delta;
		float** pW = W;
		float** pL = L;
		float* pY = Y;
		float eta_n = eta/(float)n;
		#pragma omp parallel for shared(oY, pY)
		for(int i = 0; i < n; i++) {
			oY[i] = pY[i];
		}
		int k = Permut->ReCreate(l, true);
        for(int d = 0; d < l; d++) {
			k += Permut->GetStep();
            LtoY(k);
            CompS();
			#pragma omp parallel for reduction(+: MSE) shared(eta_n, oY, pY, pL, pW)
            for(int i = 0; i < n; i++) {
                for(int j = 0; j < n; j++) {
                    delta = (float)(eta_n * ((((pL[k][i]-pY[i])* pL[k][j]) + ((pL[k][j]-pY[j])* pL[k][i]))/2.0));
                    pW[i][j] += delta;
                    MSE += delta*delta;
                }
            }
        }
        MSE = sqrt(MSE/(float)(w*l));
		#pragma omp parallel for shared(oY, pY)
        for(int i = 0; i < n; i++) {
            Y[i] = oY[i];
		}
        delete [] oY;
    }
    return MSE;
}

int THopNet::LearnMPRMSE(float MSEmin) {
    float MSE;
    int count = 0;
    RandW();
    printf("Clear W matrix.\n");
    printf("Start LearnZRP function.\n");
    printf("#Step\tMSE\t\t(d=%i\teta=%1.10f\tlearnkind=z)\n",d,eta);
    do {
        MSE = LearnMPR();
        printf("%i\t%1.16f\n",count++,MSE);
    } while (MSE > MSEmin);
    return EXIT_SUCCESS;
}

float THopNet::LearnP(void) {
    if(l == 0)
        return 0.0;
    float MSE = 0.0;
    float* nL;
    if((nL = new float[n]) == NULL) {
        /* Error */
        return EXIT_FAILURE;
        /* Error */
    } else {
        float delta = 0.0;
		int k = 0;
		float** pW = W;
		float** pL = L;
		float* pY = Y;
		float eta_n = eta/(float)n;
        for(int d = 0; d < l; d++) {
			k = d;
			int ki = k < l - 1 ? k + 1 : 0;
			LtoY(ki);
			CompS();
			#pragma omp parallel for shared(k, ki, nL, pL, pW, pY)
            for(int i = 0; i < n; i++) {
				nL[i] = 0.0;
				for(int j = 0; j < n; j++) {
					nL[i] += pW[i][j] * pY[j];
				}
				nL[i] = pL[ki][i] - nL[i];
            }
			#pragma omp parallel for reduction(+: MSE) shared(eta_n, k, nL, pL, pW)
            for(int i = 0; i < n; i++) {				
                for(int j = 0; j < n; j++) {
                    if(i != j) {
                        pW[i][j] += eta_n * nL[i] * pL[k][j];
                        MSE += (eta_n * nL[i] * pL[k][j])*(eta_n * nL[i] * pL[k][j]);
                    }
                }
            }
        }
        MSE = sqrt(MSE/(float)(w*l));
        delete [] nL;
    }
    return MSE;
}

int THopNet::LearnPMSE(float MSEmin) {
    float MSE;
    int count = 0;
    ClearW();
    printf("Clear W matrix.\n");
    printf("Start LearnP function.\n");
    printf("#Step\tMSE\t\t(d=%i\teta=%1.10f\tlearnkind=p)\n",d,eta);
    do {
        MSE = LearnP();
        printf("%i\t%1.16f\n",count++,MSE);
    } while ((MSE > MSEmin) && (count < n));
    return EXIT_SUCCESS;
}

float THopNet::LearnS(void) {
    if(l == 0)
        return 0.0;
    float MSE = 0.0;
    float* nL;
    if((nL = new float[n]) == NULL) {
        /* Error */
        return EXIT_FAILURE;
        /* Error */
    } else {
        float delta = 0.0;		
		float** pW = W;
		float** pL = L;
		float eta_n = eta/(float)n;
		int k = 0;
        for(int d = 0; d < l; d++) {
			k = d;
			#pragma omp parallel for shared(k, nL, pL, pW)
            for(int i = 0; i < n; i++) {
				nL[i] = 0.0;
				for(int j = 0; j < n; j++) {
					nL[i] += pW[i][j] * pL[k][j];
				}
				nL[i] = pL[k][i] - nL[i];
            }
			#pragma omp parallel for reduction(+: MSE) shared(eta_n, k, nL, pL, pW)
            for(int i = 0; i < n; i++) {				
                for(int j = 0; j < n; j++) {
                    if(i != j) {
                        pW[i][j] += eta_n * nL[i] * pL[k][j];
                        MSE += (eta_n * nL[i] * pL[k][j])*(eta_n * nL[i] * pL[k][j]);
                    }
                }
            }
        }
        for(int d = l - 1; d >=0; d--) {
			k = d;
			#pragma omp parallel for shared(k, nL, pL, pW)
            for(int i = 0; i < n; i++) {
				nL[i] = 0.0;
				for(int j = 0; j < n; j++) {
					nL[i] += pW[i][j] * pL[k][j];
				}
				nL[i] = pL[k][i] - nL[i];
            }
			#pragma omp parallel for reduction(+: MSE) shared(eta_n, k, nL, pL, pW)
            for(int i = 0; i < n; i++) {				
                for(int j = 0; j < n; j++) {
                    if(i != j) {
                        pW[i][j] += eta_n * nL[i] * pL[k][j];
                        MSE += (eta_n * nL[i] * pL[k][j])*(eta_n * nL[i] * pL[k][j]);
                    }
                }
            }
        }
        MSE = sqrt(MSE/(float)(2*w*l));
        delete [] nL;
    }
    return MSE;
}

int THopNet::LearnSMSE(float MSEmin) {
    float MSE;
    int count = 0;
    ClearW();
    printf("Clear W matrix.\n");
    printf("Start LearnS function.\n");
    printf("#Step\tMSE\t\t(d=%i\teta=%1.10f\tlearnkind=s)\n",d,eta);
    do {
        MSE = LearnS();
        printf("%i\t%1.16f\n",count++,MSE);
    } while ((MSE > MSEmin) && (count < n));
    return EXIT_SUCCESS;
}
